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Abstract 

Inviscid flux Jacobian matrices and their properties used in numerical solutions 
of conservation laws are extended to general, equilibrium gas laws. Exact and 
approximate generalizations of the Roe average are presented. Results are given 
for one- dimensional flow, and then extended to three-dimensional flow with time- 
varying grids. 


Introduction 

Most treatments of inviscid terms in the numerical solution of conservation laws 
utilize the properties of flux Jacobian matrices. For central difference methods, 
the Beam- Warming scheme [1] requires the true flux Jacobians matrices, and their 
eigenvalues and eigenvectors are needed for the diagonal algorithms [2]. Most up- 
wind methods, such as the Steger- Warming flux- vector splitting [3], the van Leer 
flux- vector splitting [4], and the Roe approximate Riemann solver [5], all utilize 
the properties of the flux Jacobian matrix. Their original derivations relied on the 
algebraic simplicity of the perfect gas law. Many flows of current interest involve 
departure from a perfect gas due to vibrational excitation, dissociation, and ioniza- 
tion, although the assumption of thermodynamic and chemical equilibrium is still 
valid. The purpose of this paper is to derive the flux Jacobian matrices and their 
properties for a general, equilibrium gas law, and to present the generalization of 
the Roe average used in Roe’s approximate Riemann solver. Generalizations of the 
flux- vector splitting methods are reported elsewhere [6-9]. 

The generalizations of the Roe average proposed by other investigators [10-12] 
are all based on an approximate flux Jacobian matrix, or utilize some other approx- 
imations. In the earlier stages of this work, whose results are reported in Refs. 7-9, 
it was established that a Roe-averaged state exists for an equilibrium gas, but its 
precise value is not uniquely defined. A particular method for obtaining a unique 
state was proposed. In the present paper, a new, exact definition of a unique Roe- 
averaged state is given. Since its exact implementation may not always be practical, 
several approximations to the relations are also given. For simplicity, the analysis 
is first, presented for one-dimensional flow with a fixed grid. The results are then 
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generalized to three-dimensional flow with time- varying grids. 


Flux Jacobian Matrices for One-Dimensional Flow 

The primitive variables defining a fluid state are the density p, velocity u, and the 
internal energy per unit mass t. Note that there is an arbitrary additive constant 
in the definition of €. Since conservation laws are expressed in terms of conserved 
quantities per unit volume, it is convenient to introduce the internal energy per 
unit volume t = pe . The corresponding set of conservative variables U can be 
represented by the algebraic column vector 


U = 



( 1 ) 


where m — pu is the momentum per unit volume, and e — e + \pu 2 is the total 
energy per unit volume. 

The calculation of the flux of U plays a central role in the numerical solution of 
conservation laws. The set of inviscid flux components F is given by the algebraic 
column vector 

F = 


M 


pu 

P 

— 

mu + p 

E 


eu + pu 


( 2 ) 


where M, P, and E are the flux of mass, momentum, and energy, and the pressure 
p is given by a general equation of state of the form 


P = P{P>£)- 


(3) 


The derivatives will be denoted by 

* = (t)r and ^(D; (4) 

If h = (e + p)/p is the specific enthalpy, the speed of sound c can then be expressed 
as 

c 2 = x + (5) 

It will also be convenient to define nondimensional parameters 

P<? . _ n P f C \ 

7 — and 7 = 1 + -. (o) 

V e 
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Using Eq. ( 5 ) one can also express 7 as 

7 = — — (7 + “)• ( 7 ) 

7 — K. € 

Note that while ac, c, and 7 have well defined values, the values of \ and 7 depend 
on the choice of arbitrary constant in the definition of e. 

An important special equation of state is that for a thermally perfect (but calor- 
ically imperfect) gas, which has the form 

p^pf(e). (8) 

This law is valid for a dilute gas consisting of a single chemical species, and is also 
a very good approximation for air below the temperature when oxygen starts to 
dissociate (approximately 2000 K). Using Eqs. ( 4 )-( 7 ) one can readily obtain the 
relations 

P — x{ e )p + «(e)c (9) 

and 

7 = /e + 1, 7 = 7 + -. (10) 

If one further specializes Eq. (8) by letting /(e) be just a linear function of e, one 
obtains the equation for a gas that is also calorically perfect — better known simply 
as a perfect gas. This law is valid for a gas consisting of structureless particles, and 
is also a very good approximation for air at low temperatures. The derivatives \ and 
k. are now constants. It follows from Eqs. ( 9 ) and ( 10 ) that a great simplification is 
obtained if one chooses the arbitrary constant in the definition of e so that x — 0 
and 7 = 7. In fact, this choice is tacitly made in the usual definition of a perfect 
gas. It is therefore also customary to choose the arbitrary energy constant for a 
general gas so that x approaches zero at low temperatures. 

The differential expression dF = AdU defines the flux Jacobian matrix opera- 
tor A. The differential 

dp = xdp + K>de, (11) 

can be rewritten in terms of the differentials of the conservative variables as 

dp — Ki dp — ku dm. + k de, (12) 

where K 1 = The matrix A can then be written as 

0 1 0 

A = Ah - u 2 (2 - k)u k , ( 13 ) 

(Ki—H)u H — Ku 2 (l + /c)w_ 
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where H = h + ~u 2 is the total enthalpy per unit mass. The three eigenvalues of A 
are readily found to be 

A] = u, A 2 = u + c, and A 3 = u — c. (14) 


The corresponding right, eigenvector matrix R is 


R = 



1 1 
u c u — c 
H A cu H — cu 


(15) 


where K 2 = 
form 


\u 2 ~x/k = H - c 2 /k, while the left eigenvector matrix R 1 


R “ 1 = 


r 1 -k./c 2 

\ (A'l !<? -«/c) 

.§ (A' 1 /c 2 +«/c) 


Ku/c 2 ~k/c 2 

— | (ku/c 2 — 1 /c) \k!c 2 

— | (ku/c 2 -f 1/c) \k/c 2 _ 


takes the 


(16) 


The only difference between the expressions in Eqs. (13), (15), and (16) and the 
corresponding expressions for a perfect gas is the presence of \ in the terms A 1 and 
Ko, and the fact that k, is a- variable instead of a constant. By writing A] and A 2 
in terms of H and c 2 , using Eq. (5), one can obtain expressions in which the only 
difference with the perfect gas expressions is the presence of the variable k. 

One can define functions of the matrix A through 


f{A) = /(A])P] 4- /( A 2 )P 2 A /(A 3 )P 3 , (17) 


where the projection operators P r are the tensor products 

Pi = RiRT 1 . 


(18) 


Examples of /(A) are A, |A|, sgnA = |A|/A, and A* = (A ± |A|)/2. The formula for 
Pi is 

1 —Ki/c 2 Kujc 2 —k/c 2 


2 / _2 


Pi — (l — K\/c 2 )u K,u 2 / 


■Ku/i 


Using the fact that 


[u 2 - KjH/c 2 (kH/c 2 -1)u 1 — nH / & 


Pi A P 2 A P 3 — 7, 


(19) 


(20) 


where I is the unit matrix, one can express the other two projection matrices as 


P 2 = -[±A/c A 1(1 A u/c) - Pi 
3 2 


(21) 
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One can easily establish that 


AU = F + (7 - 1 - k) P 


( 22 ) 


It follows that the necessary and sufficient condition for the homogeneity property 
F = AU is given by a thermally perfect gas (Eq. ( 8 )). Using Eq. (17), one can 
expand F for a thermally perfect gas as 


i= 1 

where each F{ is associated with A;. The final expressions for the F{ are 


A 2 p (7 - 1 ) 

1 

u 

A 2 p 

, F 2 = ^ 

1 

u ± c 

7 


3 27 

H ±cu m 


(23) 


(24) 


where 7 = 7 (e). Eqs. (23) and (24) form the basis for the generalized Steger- 
Warming flux- vector splitting [6-9]. 


Generalized Roe Average for One- Dimensional Flow 

Among the various approximate Riemann solvers, the most common one uses the 
Roe average because of its simplicity and its ability to satisfy the jump conditions 
across discontinuities exactly. The derivation in Ref. 5 for a perfect gas employed 
parameter vectors. To obtain a generalization for an equilibrium gas, a different, 
more direct approach is used here. 

In approximate Riemann solvers based on local linearization, the flux at a point 
separating two states U^and Ur is based on the eigenvalues and eigenvectors of 
some average A. The optimum choice for A is one satisfying 

AF = lAF, (25) 

where A(-) = (■)# — (-)l* This choice of A captures discontinuities exactly. One 
way of obtaining A is to seek an average state £7, such that 

A = A(U). (26) 

The notation U implies only those variables that appear explicitly in Eq. (26). Such 
a state is known as a Roe-averaged state, and was derived by Roe for a perfect gas. 
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The generalization to an equilibrium gas is obtained by substituting Eqs. (1), (2), 
(13), and (26) into Eq. (25). The second component of Eq. (25) results in 


Pr u2 r ~ Pl u2 l + Ap — 


(b- - 1 )k 2 + x 

+ K Ae + 7^ K (pR u2 R ~ PL u l )• 


A p + (2 — k)u(pru r - plul) 


(27) 


The average velocity u must be some linear combination of ur and u R . Therefore 
let 

u = aui + (3u r . (28) 

Since ul and u R can vary independently, the products u R , u 2 R , and ulur are also 
independent. After substituting Eq. (28) into Eq. (27), the terms involving each 
of those products must satisfy the equation separately. From the coefficients of v? L 
one obtains the quadratic equation 


a 2 (pR - pl) + 2 apL - PL — 0. 


(29) 


The only root that is finite when p R = pl is 

\fpi 

a = — _ 

v Pl + y/PR 


(30) 


From the coefficients of u 2 R one can show that (3 = 1 — a. The equation obtained 
from the coefficients of u^ur is also satisfied by these expressions for a and (3. 
Therefore Eq. (28) can be rewritten as 

u = ocul -f (1 — ot)uR. (31) 


Eqs. (30) and (31) are the identical relations derived by Roe for a perfect gas. The 
remaining terms in Eq. (27) result in the new condition 

% Ap -f 7c Ae = Ap. (32) 


This is just the discrete form of Eq. (11), averaged between the two states. This 
last condition is automatically satisfied for a perfect gas. In a similar manner, the 
third component of Eq. (25) results in the additional relation 

H = aH L + (1 - a)H R , (33) 

which is also true for a perfect gas. From the definition of FT, Eqs. (31) and (33) 
can be combined to define the Roe-averaged specific enthalpy as 

h = ahL + (1 - ot)h R + — a(l - a)(Au) 2 . (34) 
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Note that h could lie outside the range given by h.L and Jir if the magnitude of A u 
is sufficiently large. The Roe-averaged sound speed is given by Eq. (5) as 

c 2 = x + (35) 


For a perfect gas, Eqs. (30), (31), (33), (34), and (35) are sufficient to define uniquely 
A;, R , and R \ since x — 0 and ^ is a given constant. 

For an arbitrary equilibrium gas, Eq. (32) provides only one relation for the 
variables \ and n. We thus have the paradoxical situation that not only does a 
Roe-averaged state exist for an equilibrium gas, its precise value is not uniquely 
defined. For the special case in which state L and R are precisely those that satisfy 
the jump conditions across a discontinuity, Eqs. (30) through (35) are consistent 
with the exact Riemann solver, even though x and « are not uniquely defined. For 
a stationary shock wave, 

AM = AP - AH = 0. (36) 


By combining the definitions 
readily show that 


of M and P with Eqs. (30), (31), and (36), 


M 2 

PL PR 


ulur 


A p 
A p 


one can 
(37) 


The definitions of 714, F, and H can be combined with Eqs. (36) and (3/ ) to derive 


the relation 


At — h A p. 


(38) 


By substituting Eqs. (32) and (38) into Eq. (35), it follows that for a stationary 
shock wave, if \ and k satisfy Eq. (32), then 


c 2 


A p 
Ap 


(39) 


From Eqs. (37) and (39) it follows that for a stationary shock wave the magnitude 
of the Roe-averaged velocity is equal to the Roe-averaged sound speed. This could 
have been predicted ahead of time from Eq. (14), since one of the eigenvalues had 
to be equal to zero. From the definitions of M and H and Eq. (36) it follows that 



ApAh 
PL + PR 


Substituting Eq. (40) into Eq. (34) one obtains 


(40) 


7 - PL^L + PR^R 

h — 

PL + PR 


(1 - a) 2 h L + a 2 h R 
(1 — a) 2 + a 2 


( 41 ) 
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It is easy to show that Eqs. (39) and (41) are also valid for a non -stationary shock 
wave. The values of h and c 2 as given by Eqs. (39) and (41) will in general not be 
consistent with a thermally perfect gas law, except for the special case of a perfect 


gas. 


It is clear that unique values of x and « must be defined in terms of the ther- 
modynamic states L and R . Accurate numerical calculations for air [13] show that 
X and 7c can have a non-convex behavior if the states L and R are far apart. A 
clue to a simple definition for the two derivatives can be obtained by considering 
the special case e L = 7r, i.e., A? = 0. Integrating Eq. (11) and substituting into 
Eq. (32) yields 




x{p,e L )dp. 


(42) 


Thus x is the integrated average of x along the straight-line path between states 
L and R in the p-e plane for this special case. Actually, there is an infinite set of 
paths that can be used to define x, but the straight-line path is the simplest one 
that can be defined for an arbitrary function. Similarly, for A p = 0 one obtains 


k = 


A p 
A? 



n{pL, T) de. 


(43) 


Let an arbitrary path between any two states L and R be defined parametrically 
by the functions p{t) and ?(/), where the parameter t is normalized so that ii — 0 
and i R = 1. Integrating Eq. (11) along this path, one obtains 

A p= f x[p(i),t{i)]p'(i)di+ [ «[/>(<)» (44) 
Jo Jo 

The simplest choice is the straight-line path 

p{i) — Pl + t (45a) 

e(t) = e L + tAe. (45b) 

Substituting Eqs. (45) into Eq. (44), and comparing with Eq. (32), yields the general 
relations 


X = f x[p(t),e(i)]dt, (46a) 

Jo 

k = f K,{p(t),e(t)} dt. (46b) 

Jo 
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Eqs. (45) and (46) give unique definitions of £ and « satisfying condition (32) for 
arbitrary values of A p and A?, including the limiting case A p = Ae — 0. From 
Eqs. (34) and (35) one sees that even if A u = 0, c 2 is not equal to the integrated 
average of c 2 . Since h is a smoothly varying function (compared to x a nd k), it is 
reasonable to expect that c 2 will always be positive. For a thermally perfect gas 
law obeying Eq. (8), and Ae ^ 0, Eqs. (46a) and (46b) become 


X = PLPR Ae 



(Ae — e A p) 2 ’ 


R = PlPr Ae 



f de 

(A?- eA p) 2 ' 


(47a.) 

(47b) 


Given an equation of state, or some interpolation representation (such as Ref. 13), 
the integrals in Eqs. (46) can be evaluated for any two end states L and R. Since 
the exact evaluation may not be practical, some approximate quadratures may be 
required. Let £ and k, be approximations to Eqs. (46a) and (46b). They will not 
satisfy Eq. (32) exactly. One therefore requires values of £ and k. satisfying Eq. (32) 
which are closest to £ a ^d £. This can be formulated geometrically as projecting 
the point £, k onto the straight line defined by Eq. (32). But in order for the Roe- 
averaged state to be independent of the arbitrary constant in the definition of e, one 
must first recast the problem so that geometric relationships will not be affected by 
the choice of this constant. This can be accomplished if one first divides Eq. (32) 
by k. The slope of the straight line for the variables 1/H and £/k is now given by 
A p and A p, both of which are uniquely defined by states L and R. A further scale 
factor s with the dimension of £ must, be introduced, since £ is not dimensionless. 
The projection on to the straight line is then defined by the relation 


A p + zz A p — — A p -f — A p. 

hi R /C K 

(48) 

If one introduces the error 


6p — A p — x A p — k Ae 

(49) 

and the quantity 

D = (sAp) 2 + (A pf, 

(50) 

one can solve Eqs. (32) and (48) to obtain the final relations 


_ Dx + ^ 2 ApSp 
X ~ D-ApSp ’ 

(51a) 

_ Dk 

D — A pSp 

(51b) 
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A natural choice for the scale factor s is 

5 — c 2 = x + K'h, (52) 

where the same quadrature approximation that was used to calculate £ and £ is 
assumed. Note that the s given by Eq. (52) is guaranteed to be positive. This is 
not necessarily true if one defined s' = £ + 

Let pm - {pL + Pr)I 2 and e M = (e L + 7r)/ 2 define the midpoint state M. Then 
possible quadrature rules for £ are the rmdpoint rule 

X~Xm, (53) 

the trapezoidal rule 

X = (XL + Xi?)/ 2 > ( 54 ) 

and Simpson’s rule 

X ~ (XL + 4 Xm + Xfl)/ 6 ? (55) 

with analogous formulas for 7c and c 2 . Eqs. (53) and (54) are exact if p is a quadratic 
polynomial in p and ?, while Eq. (55) is exact if p is a cubic polynomial. While the 
expressions for £ and 7c given in Refs. 7-9 used approximation (54), they differ from 

the present results since Eq. (48) was written in terms of 1 /ac and x//c. If the states 
L and R are reasonably close, approximations (53) or (54) should be adequate. 
For large separation of the two states, Simpson’s rule (55) may be required. The 
quantity D/p 2 L is a useful nondimen si onal parameter measuring the separation of 
the two states. 

While the above relations are all that are required to construct a Riemann solver 
using Roe’s linearization, an additional algebraic simplicity can be achieved by 
expressing differences in conservative variables in terms of differences in primitive 
variables. If one formally defines 


P = y/PLpR, 

(56) 

one obtains the identities 


A (pu) — pAu -f uAp 

(57) 

and 


A (pu 2 ) — 2puAu -fi u 2 Ap. 

(58) 


An important quantity in the approximate Riemann solver is the column vector 
R~ 1 AU. Its components are the jumps in the characteristic variables. Using 
Eqs. (32), (57), and (58), it can be expressed simply in term of A p, A p and A u as 

A p — Ap/c 2 

R^AU = \{Aplc 2 + p Au/c) . (59) 

_ \(Ap/c 2 — p Au/c) _ 
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These expressions have the same form as for the perfect gas case. The quantities x 
and k do not appear explicitly. 


Comparisons with Other Formulations 

The formulas for Roe’s approximate Riemann solver described by Grossman and 
Walters [10] involve three approximations. They assume that 1 ) A 7 and A 7 are 
small, 2 ) A 7 can be related to A p via an isentropic formula, and 3) 7 is nearly 
equal to 7 . All these approximations can be poor if the two states are far apart. If 
one notes that 

c 2 = Az 21 7 h, (60) 

7 

their final expressions are equivalent to setting x — 0 ? with 7 and 7 replaced by 
their arithmetic averages. The calculations of Ref. 13 show that setting x = 0 in 
the first eigenvector of Eq. (15) may not be justified. 

Glaister [ 11 ] follows an analysis similar to the present one, except that he employs 
e instead of e in Eqs. (3), (4), and (11). As a result it is necessary to introduce p 
via Eq. (56), as well as to define 

e - ae L + (1 - a)e R (61) 

in order to obtain Eq. (32) in terms of Ae. In the limit of a perfect gas, note the 
inconsistency between Eqs. (34) and (61). He suggests a different integration path 
in order to define his x and ft. Let Pl^r define the fictitious state T, and pR,e L 
define state B. Then x and ^ are defined as the arithmetic averages of the values 
resulting from integrating Eq. (32) via paths L-T-R and L-B-R , respectively. This 
results in 


x = ^(Ap + PB - Pt)/Ap, 

(62a) 

k = ^(A p + PT - pb)/ Ae. 

(62b) 


These are replaced by x = \{XL + Xr) if Ap = 0, and 7c — \(kl + kr) if Ae = 0. 
If the states L and R are far apart, and the equation of state is non-convex, the 
introduction of the fictitious states T and B could give poor results. 

Liou et al [12] also follow Glaister’s analysis, employing e instead of e, but suggest 
using the state M defined by Eqs. (56) and (61) instead of M to calculate their 
approximate x an( i « by the midpoint rule (53). From the interpretation of x an d 
k as quadrature approximations to integrals, it would appear that evaluating at the 
average state M would give a better approximation in general. Their formulas for 
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the values of x and * satisfying Eq. (32) which are closest to \ an d * a * so break 
down when either A p or Ae approaches zero. 


Generalization to Three-Dimensional Flow 

The generalization of the above results to three-dimensional flow can be accom- 
plished in a compact manner by employing the vector approach of the author [14]. 
Here the word vector refers to a physical vector such as velocity or momentum, as 
distinguished from an algebraic, vector representing a set of variables. If u is the 
fluid velocity vector, then U can be represented by the column vector 



(63) 


where m = pu is the momentum vector per unit volume, and e — e + |pu • u is the 
total energy per unit volume. Let n be the unit normal vector in a positive direction 
to a cell surface in a finite- volume grid, or a coordinate surface in a finite-difference 
grid. If v n is the normal component of the velocity of a time-varying surface, and 
u n = u • n, one can define the normal relative velocity component u' = u n - v n . 
The set of inviscid normal flux components F n is given by the column vector 



~ M~ 


i 

i 

F n = 

P 

— 

mid + pn 


E 


_ eu 1 A pu n _ 


(64) 


where M, P, and E are the normal flux of mass, momentum, and energy. 

The flux Jacobian requires expressing dF n in terms of dJJ . For a fixed v n and n, 
the first component of dF n can be written as 


d(pu') — —v n dp + n • dm. 

This can be rewritten in the form of a matrix multiplication as 


(65) 


d(pu') — [ —v n 


n 


0 


dp 

dm 

de 


( 66 ) 


where the dot product is implied in multiplying the second element of the row 
vector by the second element of the column vector. Applying the same procedure 
to the other components of dF n , one can define the flux Jacobian matrix operator 
A satisfying dF n = AdU , using the convention that in forming the product of a 
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matrix element with a vector element, a dot product is implied if each element is 
either a physical vector (e.g. u or n) or a tensor (e.g. un or nu). The matrix A 
can then be written as 

-v n n 0 

A = Kjn — u n u un — /am + u'l /cn , (67) 

(K 1 —H)u n Hu-ku^u u' + nu n _ 

where AT = * u + \-> H = h + |u • u is the total enthalpy per unit mass, and I 

is the identity tensor. The three eigenvalues of A are 

Aj = u , Ao = u + c, and A3 = u — c. (68) 

Aj is a repeated eigenvalue, requiring a set of linearly independent eigenvectors, 
which can be defined with full generality in term of an arbitrary set of spatial basis 
vectors a;, and a set of reciprocal basis vectors a- 7 satisfying a 2 - a- 7 = where 8j 
is the Kronecker delta. One can then define a n2 = n • a*, b 2 = n x a 2 , a J n — n • a- 7 , 
and b J = n x a- 7 . If (3 is an arbitrary scalar, the right eigenvector matrix R can be 
written in the most general form as 

Oni I 1 

R = a ni u + f3bj u + cn u-cn , (69) 

a m AT + (3h i • U H + CU n H — CU-n 

where AT = ^u • u - ^//c, while the left eigenvector matrix P~ 3 takes the form 

' a{ (l - AT/c 2 ) - b J • \i/j3 a J n Ku/c 2 +b j //3 -a J n n/c 2 ' 

R~ l = \ (AT/c 2 - u n /c) (/cu/c 2 - n/c) \k/c 2 . (70) 

\ (AT/c 2 + u n /c) (ku /c 2 + n/c) |/c/c 2 _ 

Note that i and j take on values from 1 to 3, so that R has five columns and A -1 
has five rows. A useful choice for the basis vectors is to let one of the a; be parallel 
to n, so that the corresponding b; = 0. The remaining a; are then chosen to lie in 
the plane perpendicular to n, so that their corresponding a ni = 0. 

Eq. (17) which defines a function of A is still valid, with P a given by 

l-A]/c 2 /cu/c 2 — /c/c 2 

Pj = u n n-Aiu/c 2 kuu/ c 2 — nn -j- 1 —/cu/c 2 , (71) 

u n 2 — KiH/c 2 kHu/c 2 — u n n 1 — /c H/c 2 _ 

and the other two projection matrices given by 

P, = 1[±A/c + /(1t«7c)-J>,]. (72) 

3 L 
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In the diagonal algorithms [2] one encounters the product of two eigenvector 
matrices, each of which is associated with a different surface. Let the surface with 
the unit normal vector rife have associated with it the set of basis vectors a i k with 
reciprocal basis vectors aj, , and an arbitrary scalar ft. One can then define a n ;*. = 
n k * a ifc , b lk = n*. x a ik , a 3 nk = n k • a.{ , and b[ = n k x aj. Similarly, let the surface 
with the unit normal vector n; have associated with it the set of basis vectors a^ 
with reciprocal basis vectors a j , and an arbitrary scalar ft. One can then define 
anil = n I • Siil, b u - rn x a ih a J nl = ni-a{, and bj = m x aj. Then the five by five 
matrix Ri can be written as 


R k J Ri 


a nil a i k + b u • b ; fc ft/ft 

|ftb« • n k /c 

-fftbii • njt/c 


cb£ • n z/ft 
|(1 + • n/) 

|(1 - n k • nz) 


-cbj. • n i/ft. 
|(1 - • n z ) 

|(1 + n*. • ii;) _ 


(73) 


From the form of (73) it follows that the matrix only depends on geometric factors 
if both ft and ft are chosen to be proportional to the sound speed c. 

For a thermally perfect gas one can expand F n as 


F — F 

■L n — / 1 n 1 ■> 


(74) 


l~ 1 


where the F ni are given by 
ft p(l - 1) 


F n l ~ 


7 


1 

u 


L 2-u • u f - TfTT) j 


F 2 = 


A 2 p 

3 

27 


1 

u ± cn 

H ± cu n 


(75) 


The generalization of the Roe average for a surface separating Ul and Ur is 
obtained by substituting Eqs. (63), (64), (67), and (26) into A F n = AAU , with n 
and v n fixed at the surface. In the momentum equation, u must be some linear 
combination of the vectors uj,, u.r, and n. Since these three vectors can vary 
independently, the dot products Ul • u^, u L • u r, u r • Ur, uj, • n, and u# • n are 
also independent. In a manner similar to the one outlined for the one-dimensional 
case, one can easily show that 

u = aui + (1 - a)ufl, (76) 

where a is still given by Eq. (30). Eqs. (32), (33), and (35) are still valid, with h 
given by 

h = ahi + (1 — ct)h.R -f -a(l — a)Au • Au. (77) 

Z 
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The previously derived expressions for x and K remain unchanged. The expression 
for R~ 2 AU becomes 


R~ J AU = 


’ a J n (Ap - Ap/c 2 ) + p(b 3 • Au) /j3 
f (A p/c 2 + pn • Au/c) 

|(A p/c 2 - pn • Au/c) 


(78) 


Conclusions 

Inviscid flux Jacobian matrices and their properties which are useful in CFD 
techniques have been derived for an arbitrary, equilibrium gas. Both exact and 
approximate generalizations of the Roe average have been presented. 
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